# estimate John Jackson's model using multinomRob
#

options(width=130);

library(multinomRob)

# read John Jackson's replication dataset (PA 2002, Winter)
#
jjdat <- read.table(file="pareplic.asc", header=T);

#> names(jjdat)
# [1] "Voivodship"    "ID"            "UD.KLD"        "SLD"
# [5] "PSL"           "UP"            "Catholic"      "Other"
# [9] "Constant"      "New.Jobs"      "Unemp"         "SOE"
#[13] "Church.Att"    "Schooling"     "Age"           "Farmer"
#[17] "Village"       "Suchoka.UD"    "Pawlak.PSL"    "Kwasn.SLD"
#[21] "Zioqkowska.UP" "Voters"

jjuse <- jjdat;
#total number of voters
jjuse$TotalVotes   <- jjdat$Voters * 1000;
jjuse$SLD      <- jjdat$SLD/100      *jjuse$TotalVotes;
jjuse$PSL      <- jjdat$PSL/100      *jjuse$TotalVotes;
jjuse$UP       <- jjdat$UP/100       *jjuse$TotalVotes;
jjuse$Catholic <- jjdat$Catholic/100 *jjuse$TotalVotes;
jjuse$Other    <- jjdat$Other/100    *jjuse$TotalVotes;
jjuse$UD.KLD   <- jjdat$UD.KLD/100   *jjuse$TotalVotes;

jjuse$Age  <-  jjdat$Age /10;

## Set parameters for Genoud
zz.genoud.parms <- list( pop.size             = 2000,
                        wait.generations      = 30,
                        P5                    = 500,
                        scale.domains         = 2
                        )

# estimate the model as done in the Mebane-Sekhon AJPS 2004 paper

load("jj.RData");
starts <- muljj$genoud$par;
#starts[c(8,18,26,41)] <- muljj$mtanh$coeffvec[c(8,18,26,41)];

muljj <-
  multinomRob(
    list(SLD ~      New.Jobs + Unemp + SOE + Church.Att + Schooling + Age + Kwasn.SLD ,
         PSL ~      New.Jobs + Unemp + SOE + Church.Att + Schooling + Age + Farmer + Village + Pawlak.PSL ,
         UP ~       New.Jobs + Unemp + SOE + Church.Att + Schooling + Age + Zioqkowska.UP ,
         Catholic ~ New.Jobs + Unemp + SOE + Church.Att + Schooling + Age ,
         Other ~    New.Jobs + Unemp + SOE + Church.Att + Schooling + Age ,
         UD.KLD ~   0 + Suchoka.UD ),
    jjuse,
    starting.values = starts,
    genoud.parms = zz.genoud.parms,
    print.level = 3);

summary(muljj, weights=TRUE);

save(muljj, file="jj.a.RData")
